Example data analysis for manuscript

baRulho:quantifying habitat-induced degradation of (animal) acoustic signals

Author

Marcelo Araya-Salas

Published

August 10, 2023

0.1 Load package

Code 
library(baRulho)
library(warbleR)

0.2 Synthetize sounds

Create synthetic sounds to be used for making the master sound file for playback experiments:

Code 
synth_data <-
    synth_sounds(
        replicates = 3, # number of replicates for each unique combination of varying features
        frequencies = seq(0.5, 10, length.out = 20),
        durations = c(0.2, 0.1),
        am = TRUE, # amplitude modulation
        fm = TRUE, # frequency modulation
        sig2 = 0.8, # frequency modulation parameter
        shuffle = TRUE # randomize the position of sounds 
    )

The output is of class data frame and extended selection table (warbleR package format, here printed as a data frame):

Code 
head(synth_data, 10)
[1] 68 70 72 119 78 120 86 235 33 93 140 16 80 229 48 375 192 370 [19] 44 9 67 58 127 276 168 227 207 355 225 107 275 406 242 343 224 407 [37] 249 237 387 263 175 117 381 377 295 322 273 368 203 76 122 43 92 52 [55] 55 121 18 277 130 158 139 83 145 155 294 313 202 285 136 126 300 82 [73] 291 233 310 177 111 435 421 328 278 248 239 282 352 408 386 426 213 74 [91] 374 231 409 422 210 240 412 106 380 437 360 357 267 415 262 372 396 379 [109] 400 129 113 64 144 14 99 79 71 236 108 8 6 148 19 73 165 167 [127] 221 264 115 38 178 3 271 284 230 228 132 371 404 369 269 299 162 197 [145] 365 376 425 329 255 152 414 307 402 319 173 303 347 91 429 321 391 327 [163] 287 397 247 433 385 317 417 84 60 45 100 34 134 63 193 125 23 15 [181] 94 53 204 154 344 281 353 138 42 56 85 211 289 201 217 309 434 241 [199] 137 256 392 220 411 214 419 182 358 333 292 250 59 364 174 351 216 362 [217] 293 131 1 27 69 146 7 35 11 194 123 10 65 114 21 36 222 169 [235] 186 279 75 96 180 98 232 286 195 170 135 346 373 324 160 226 254 252 [253] 337 345 367 331 166 95 416 251 315 366 301 270 389 12 403 323 320 427 [271] 290 410 171 388 325 390 420 142 141 62 47 128 29 26 133 185 147 13 [289] 31 118 49 159 190 172 219 302 104 81 208 25 260 297 314 288 149 438 [307] 394 418 296 206 184 238 363 336 354 356 188 57 424 215 423 341 283 274 [325] 378 28 405 326 339 413 304 428 187 361 338 335 430 54 112 150 151 212 [343] 110 109 359 30 306 40 156 89 77 39 102 189 305 124 50 340 234 245 [361] 266 258 198 311 268 348 103 432 265 157 431 399 37 66 196 259 436 223 [379] 398 199 46 349 209 401 280 312 205 32 90 24 105 191 5 22 342 101 [397] 261 87 97 88 143 2 41 20 183 61 4 179 298 51 164 200 161 350 [415] 246 316 153 395 181 257 330 334 318 17 244 176 308 253 116 332 393 383 [433] 243 163 272 384 218 382
sound.files selec start end bottom.freq top.freq frequency duration frequency.modulation amplitude.modulation treatment replicate sound.id
synthetic_sound_68 1 0.05 0.25 8.896 9.584 9 0.2 fm am dur:0.2;freq:9;fm;am 1 dur:0.2;freq:9;fm;am_1
synthetic_sound_70 1 0.05 0.25 3.215 4.219 4 0.2 fm am dur:0.2;freq:4;fm;am 1 dur:0.2;freq:4;fm;am_1
synthetic_sound_72 1 0.05 0.25 0.244 1.399 1 0.2 fm am dur:0.2;freq:1;fm;am 1 dur:0.2;freq:1;fm;am_1
synthetic_sound_119 1 0.05 0.25 6.209 7.184 7 0.2 fm am dur:0.2;freq:7;fm;am 1 dur:0.2;freq:7;fm;am_1
synthetic_sound_78 1 0.05 0.25 9.007 9.777 9.5 0.2 fm am dur:0.2;freq:9.5;fm;am 1 dur:0.2;freq:9.5;fm;am_1
synthetic_sound_120 1 0.05 0.25 3.918 4.800 4.5 0.2 fm am dur:0.2;freq:4.5;fm;am 1 dur:0.2;freq:4.5;fm;am_1
synthetic_sound_86 1 0.05 0.25 2.899 3.773 3 0.2 fm am dur:0.2;freq:3;fm;am 1 dur:0.2;freq:3;fm;am_1
synthetic_sound_78 2 0.05 0.25 9.007 9.777 9.5 0.2 fm am dur:0.2;freq:9.5;fm;am 2 dur:0.2;freq:9.5;fm;am_2
synthetic_sound_33 1 0.05 0.25 0.954 2.230 2 0.2 fm am dur:0.2;freq:2;fm;am 1 dur:0.2;freq:2;fm;am_1
synthetic_sound_93 1 0.05 0.25 6.006 6.790 6.5 0.2 fm am dur:0.2;freq:6.5;fm;am 1 dur:0.2;freq:6.5;fm;am_1

1 Create master sound file

This step puts all sounds together into a single sound file:

Code 
master_annotations <- master_sound_file(X = synth_data, # synthetic sound data
                  file.name = "master", # name of the sound file
                  gap.duration = 0.2) # duration of silence in between sounds

The output file is saved in the current working directory (can be modified using argument ‘path’). A similar file was used for the playback experiments detailed in the paper. The following section shows how to access the test (re-recorded) files.

These are the annotations for the sounds in the master sound files:

Code 
head(master_annotations)
sound.files selec start end bottom.freq top.freq sound.id
master.wav 1 1.000000 2.000000 1.333333 2.666667 start_marker
master.wav 2 2.050000 2.250023 7.875000 8.805000 dur:0.2;freq:9;fm;am_1
master.wav 3 2.300023 2.500045 3.208000 4.069000 dur:0.2;freq:4;fm;am_1
master.wav 4 2.550045 2.750068 0.422000 1.223000 dur:0.2;freq:1;fm;am_1
master.wav 5 2.800068 3.000091 6.905000 7.917000 dur:0.2;freq:7;fm;am_1
master.wav 6 3.050091 3.250113 8.417000 9.416000 dur:0.2;freq:9.5;fm;am_1

1.1 Download data

This code downloads the test files. The files were re-recorded during a transmission experiment at 10, 30, 65 and 100 m:

Code 
path_to_files <- "PATH_TO_FILES"  # add folder path to keep master and test files

# directory path where supplementary files have been saved

options(sound.files.path = path_to_files)

download.file("https://figshare.com/ndownloader/files/41905809", destfile = file.path(path_to_files,
    "degrad_exp_files.zip"), method = "wget")

unzip(file.path(path_to_files, "degrad_exp_files.zip"), exdir = file.path(path_to_files))

2 Find markers

The code below finds the position of the start and end markers in the test files:

Code 
# directory path where supplementary files have been saved
options(sound.files.path = "PATH_TO_FILES")

markers_in_tests <- find_markers(X = master_annotations)  # annotations of sounds in master file

head(markers_in_tests)
sound.files selec start end scores marker time.mismatch
trnsc1_100m_closed.wav 1 6.049203 7.049203 0.2227792 start_marker 0.0300544
trnsc1_100m_closed.wav 2 257.700686 258.700686 0.2726756 end_marker NA
trnsc1_100m_open.wav 3 29.502300 30.502300 0.7150524 start_marker 0.0131017
trnsc1_100m_open.wav 4 281.136830 282.136830 0.6083450 end_marker NA
trnsc1_10m_closed.wav 5 39.452253 40.452253 0.7046298 start_marker 0.0114997
trnsc1_10m_closed.wav 6 291.085182 292.085182 0.8364759 end_marker NA

The column ‘time.mismatch’ compares the time difference between the two templates on test-files against that in the master sound file. In a perfect marker detection the value must be 0, meaning that the time in between markers in the master and test files is exactly the same. In this case the average mismatch is of 14 ms and the highest of 32 ms:

Code 
# average mismatch
mean(markers_in_tests$time.mismatch, na.rm = TRUE)
[1] 0.01438455
Code 
# maximum mismatch
max(markers_in_tests$time.mismatch, na.rm = TRUE)
[1] 0.03171009

Modifing detection parameters as spectrogram type (‘type’ argument), time window overlap (‘ovlp’ argument) and hop size (‘hop.size’ argument) can be adjusted in order to improve precission. Note that for aligning all other sounds only the marker with the highest correlation will be used. Therefore the time mismatch is likely to be lower in the aligned test sounds.

3 Align sounds

Once we know the position of markers we can compute the position for all other sounds in the test files (i.e. align):

Code 
aligned_tests <-
    align_test_files(X = master_annotations, # annotations of sounds in master file
                     Y = markers_in_tests, # position of markers in test files
                     output = "data.frame") 


head(aligned_tests)
sound.files selec start end bottom.freq top.freq sound.id marker
trnsc1_100m_closed.wav 1 7.129257 7.329280 7.875 8.805 dur:0.2;freq:9;fm;am_1 end_marker
trnsc1_100m_closed.wav 2 7.379280 7.579303 3.208 4.069 dur:0.2;freq:4;fm;am_1 end_marker
trnsc1_100m_closed.wav 3 7.629303 7.829325 0.422 1.223 dur:0.2;freq:1;fm;am_1 end_marker
trnsc1_100m_closed.wav 4 7.879325 8.079348 6.905 7.917 dur:0.2;freq:7;fm;am_1 end_marker
trnsc1_100m_closed.wav 5 8.129348 8.329371 8.417 9.416 dur:0.2;freq:9.5;fm;am_1 end_marker
trnsc1_100m_closed.wav 6 8.379371 8.579393 4.171 4.839 dur:0.2;freq:4.5;fm;am_1 end_marker

This plot shows the aligned re-recorded sounds for transcect 1:

Code 
graphics.off()

lf <- rep(c(0.06, 0.5), each = 4)
rg <- rep(c(0.5, 0.95), each = 4)
horiz <- seq(0.95, 0.075, length.out = 5)
btm <- rep(horiz[-1], 2)
tp <- rep(horiz[-length(horiz)], 2)

m <- cbind(lf, rg, btm, tp)

lf <- c(rep(0.95, each = 4), 0.06, 0.5, 0, 0)
rg <- c(rep(1, each = 4), 0.5, 0.95, 0.05, 1)
horiz <- seq(0.95, 0.075, length.out = 5)
btm <- c(horiz[-1], 0.95, 0.95, 0.075, 0)
tp <- c(horiz[-length(horiz)], 1, 1, 0.95, 0.075)

m2 <- cbind(lf, rg, btm, tp)
m <- rbind(m, m2)


for (e in paste0("trnsc", 1:3)) {

    png(filename = paste0("./output/spectrograms_by_habitat_and_distance_",
        e, ".png"), res = 300, width = 4000, height = 3000)

    ss <- split.screen(figs = m)

    # # testing layout screens for(i in 1:nrow(m)) {screen(i)
    # par( mar = rep(0, 4)) plot(0.5, xlim = c(0,1), ylim =
    # c(0,1), type = 'n', axes = FALSE, xlab = '', ylab = '',
    # xaxt = 'n', yaxt = 'n') box() text(x = 0.5, y = 0.5,
    # labels = i) } close.screen(all.screens = T)

    ovlp <- 1
    colls <- seq(-110, 0, 5)
    wl <- 512

    lab_bg <- viridis(10, alpha = 0.25)[8]

    lab_bg <- "#3A3B39"
    files <- c("_10m_open.wav", "_30m_open.wav", "_65m_open.wav",
        "_100m_open.wav", "_10m_closed.wav", "_30m_closed.wav", "_65m_closed.wav",
        "_100m_closed.wav")

    files <- paste0(e, files)

    titles <- gsub(paste0(e, "_|_|.wav"), " ", files)

    # frequency label
    screen(15)
    par(mar = c(0, 0, 0, 0), new = TRUE)
    plot(1, frame.plot = FALSE, type = "n", yaxt = "n", xaxt = "n")
    text(x = 0.9, y = 1, "Frequency (kHz)", srt = 90, cex = 1.6)

    # time label
    screen(16)
    par(mar = c(0, 0, 0, 0), new = TRUE)
    plot(1, frame.plot = FALSE, type = "n", yaxt = "n", xaxt = "n")
    text(x = 1, y = 0.75, "Time (s)", cex = 1.6)

    for (i in seq_along(files)) {
        # print(i)
        screen(i)
        par(mar = c(0, 0, 0, 0))

        warbleR:::spectro_wrblr_int2(wave = readWave(file.path(path_to_files,
            files[i]), from = min(aligned_tests$start[aligned_tests$sound.files ==
            files[i]]) - 1.2, to = min(aligned_tests$start[aligned_tests$sound.files ==
            files[i]]) + 4, units = "seconds"), collevels = colls,
            ovlp = ovlp, wl = wl, flim = c(0.1, 10.6), palette = viridis,
            axisX = FALSE, axisY = FALSE, grid = FALSE)

        # add frequency axis
        if (grepl("open", aligned_tests$sound.files[aligned_tests$sound.files ==
            files[i]]))
            axis(2, at = c(seq(2, 10, 2)))

        # add time axis
        if (grepl("100m", aligned_tests$sound.files[aligned_tests$sound.files ==
            files[i]]))
            axis(1)

        lns <- c(aligned_tests$start[aligned_tests$sound.files ==
            files[i]] - min(aligned_tests$start[aligned_tests$sound.files ==
            files[i]]), aligned_tests$end[aligned_tests$sound.files ==
            files[i]] - min(aligned_tests$start[aligned_tests$sound.files ==
            files[i]])) + 1.2

        lns <- c(lns, min(lns) - 0.1, min(lns) - 1.05)

        abline(v = lns, col = "white", lty = 3, lwd = 1.2)
        abline(v = lns, col = "white", lty = 3, lwd = 1.2)
    }

    vlabs <- paste(c(10, 30, 65, 100), "m")

    par(mar = c(0, 0, 0, 0), bg = lab_bg, new = TRUE)
    # add vertical labels
    for (i in 9:12) {
        screen(i)
        # par(mar = c(0, 0, 0, 0))
        par(mar = c(0, 0, 0, 0), bg = lab_bg, new = TRUE)
        plot(1, frame.plot = FALSE, type = "n", yaxt = "n", xaxt = "n")
        text(x = 1, y = 1, vlabs[i - 8], srt = 270, cex = 1.6, col = "white",
            font = 2)
        box()
    }

    hlabs <- c("Open habitat", "Closed habitat")
    for (i in 13:14) {
        screen(i)
        par(mar = c(0, 0, 0, 0), bg = lab_bg, new = TRUE)
        plot(1, frame.plot = FALSE, type = "n", yaxt = "n", xaxt = "n")
        text(x = 1, y = 1, hlabs[i - 12], font = 2, cex = 1.6, col = "white")
        box()
    }

    dev.off()
}

Spectrograms

4 Measuring degradation

Must degradation metrics involve comparing tests sounds that were recorded at different distances from the speaker, to their reference, which is typically recorded at 1m. Hence, a column indicating the distance at which each sound was recorded is needed. In this case the recording distance can be extracted from the sound file name:

Code 
# get distance
aligned_tests$distance <- sapply(strsplit(aligned_tests$sound.files,
    "_"), "[[", 2)

# make it a numeric column
aligned_tests$distance <- as.numeric(gsub("m", "", aligned_tests$distance))

Once the distance is included in the annotations degradation metrics can be obtained with few lines of code. For instance the following code computes excess attenuation, signal-to-noise ratio, blur ratio and tail-to-signal ratio:

Code 
# run as a pipe
aligned_tests <- aligned_tests |>
    excess_attenuation() |>
    signal_to_noise_ratio(mar = 0.1) |>
    blur_ratio() |>
    tail_to_signal_ratio(mar = 0.1)  # mar = margin where to measure tail or noise

head(aligned_tests)
[1] “assuming all sound files have the same sampling rate” [1] “assuming all sound files have the same sampling rate” Preparing data for analysis (step 1 out of 3): Calculating amplitude envelopes (step 2 out of 3): Calculating blur ratio (step 3 out of 3): [1] “assuming all sound files have the same sampling rate”
sound.files selec start end bottom.freq top.freq sound.id marker distance reference excess.attenuation signal.to.noise.ratio blur.ratio tail.to.signal.ratio
trnsc1_100m_closed.wav 1 7.129257 7.32928 7.875 8.805 dur:0.2;freq:9;fm;am_1 end_marker 100 trnsc1_1m_open.wav-1 -10.934431 0.3619688 0.2684761 -0.705663
trnsc1_100m_open.wav 1 30.552300 30.75232 7.875 8.805 dur:0.2;freq:9;fm;am_1 start_marker 100 trnsc1_1m_open.wav-1 16.052347 0.7887214 0.2148822 -1.235814
trnsc1_10m_closed.wav 1 40.513753 40.71378 7.875 8.805 dur:0.2;freq:9;fm;am_1 end_marker 10 trnsc1_1m_open.wav-1 4.166308 22.8338669 0.1168380 -25.819426
trnsc1_10m_open.wav 1 31.154559 31.35458 7.875 8.805 dur:0.2;freq:9;fm;am_1 end_marker 10 trnsc1_1m_open.wav-1 6.759965 23.6130624 0.1643497 -33.034854
trnsc1_1m_open.wav 1 101.175861 101.37588 7.875 8.805 dur:0.2;freq:9;fm;am_1 end_marker 1 trnsc1_1m_open.wav-1 NA 26.9087306 NA -43.379831
trnsc1_30m_closed.wav 1 31.784472 31.98449 7.875 8.805 dur:0.2;freq:9;fm;am_1 end_marker 30 trnsc1_1m_open.wav-1 5.981216 2.5834933 0.1646304 -3.672095

Session information

R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] viridis_0.6.3      viridisLite_0.4.2  tidyr_1.1.3        baRulho_1.1.0     
 [5] ohun_1.0.0         warbleR_1.1.29     NatureSounds_1.0.4 seewave_2.2.1     
 [9] tuneR_1.4.4        rprojroot_2.0.3    formatR_1.11       knitr_1.43        
[13] kableExtra_1.3.4   remotes_2.4.2     

loaded via a namespace (and not attached):
 [1] httr_1.4.4         jsonlite_1.8.7     brio_1.1.3         assertthat_0.2.1  
 [5] yaml_2.3.7         pillar_1.9.0       backports_1.4.1    glue_1.6.2        
 [9] digest_0.6.33      checkmate_2.2.0    rvest_1.0.3        Sim.DiffProc_4.8  
[13] colorspace_2.1-0   htmltools_0.5.5    pkgconfig_2.0.3    purrr_1.0.0       
[17] scales_1.2.1       vdiffr_1.0.5       webshot_0.5.4      svglite_2.1.0     
[21] dtw_1.23-1         tibble_3.2.1       proxy_0.4-27       generics_0.1.3    
[25] ggplot2_3.4.2      pbapply_1.7-2      cli_3.6.1          magrittr_2.0.3    
[29] evaluate_0.21      fansi_1.0.4        MASS_7.3-60        xml2_1.3.5        
[33] class_7.3-22       tools_4.1.0        lifecycle_1.0.3    stringr_1.5.0     
[37] fftw_1.0-7         munsell_0.5.0      compiler_4.1.0     Deriv_4.1.3       
[41] e1071_1.7-13       signal_0.7-7       systemfonts_1.0.4  rlang_1.1.1       
[45] classInt_0.4-9     units_0.8-2        grid_4.1.0         RCurl_1.98-1.12   
[49] rstudioapi_0.14    rjson_0.2.21       htmlwidgets_1.5.4  igraph_1.5.0.1    
[53] bitops_1.0-7       rmarkdown_2.23     testthat_3.1.10    gtable_0.3.3      
[57] DBI_1.1.3          R6_2.5.1           gridExtra_2.3      dplyr_1.0.10      
[61] fastmap_1.1.1      utf8_1.2.3         KernSmooth_2.23-21 stringi_1.7.12    
[65] parallel_4.1.0     Rcpp_1.0.11        vctrs_0.6.3        sf_1.0-14         
[69] png_0.1-7          tidyselect_1.2.0   xfun_0.40